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A great deal of digital topographic information and gravity 
data exist that cover the continental United States. For ^instance, 
data sets that can be purchased from the National Oceanic and 
Atmospheric Adminstration (NOAA) include: (a) elevations 'averaged 

over 30 seconds of latitude and longitude for the cnterminus United 
States and (b) approximately a half million gravity readings, where 
the stations at*e distributed across the entire country, although in 
an irregular manner. Both the topography and the gravity data sets 
are fundamental for understanding the structure of the crust and 
lithosphere, and especially the relationships between topography and 
structure (McGinnis et al., 1979; McNutt, 1981; and others). 

The intent of this work is to show the applicability of fairly 
standard image processing techniques to processing arid analyzing 
large geologic data sets. In particular, we have utilized image 
filtering techniques to interpolate between gravity station 
locations to produce a regularly spaced data array that preserves 
detail in areas with good coverage, and that produces a continuous 
tone image rather than a contour map. We have used standard image 
processing techniques to digitally register and overlay topographic 
and gravity data, and we have displayed the data in ways that 
emphasize subtle but pervasive structural features. In this paper 
we discuss the techniques used, briefly describe the products, and 
Illustrate the potential of the methods by discussing subtle linear 
structures that appear in the processed data between the 
midcontinent gravity high and the Appalachians. 
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DATA PROCESSING TECHNIQUES 


Data processing was conducted on a PDP-11/34 minicomputer with 
interactive Image storage and display peripherals. Software Is 
based on mini-VICAR, with a variety of applications program to 
conduct filter1n(j, enhancement, and geometric operations. The 
reader Is referred to such standard references as Moik (1980) for 
further Information regarding general techniques used In Image 
processing. Our software was developed partly at the Jet Propulsion 
Laboratory, partly at the U.S. Geological Survey, Flagstaff, 

Arizona, and partly by us. 

For the purposes of this paper, the dynamic range of both the 
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range occupied by a byte variable (I.e., 8 bits or 256 discrete 
values). The byte-encoded data were then stored In data arrays, 
where the array element spacing was 30 seconds in both latitude and 
longitude. For the topography, this procedure produced a regular 
array of byte-encoded data that had 2830 rows (latitude) and 7080 
columns (longitude). The data arrays covered a latitude range from 
25* to 49* north and a longitude range from 66* to 125* west. 

The resulting topography data array could be displayed as an 
image If the byte-encoded data were converted to color or brightness 
values. An alternative and quite useful method of depicting 
topography is to generate a digital shaded relief map. Figure 1 
shows such a product that was generated on a version of the data 
file that was. first transformed to a Mercator projection, digitally 
reduced in size by a factor of seven, and then used to generate a 
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relief map. The shading process 1s based or assuming a photometric 
function for the surface (Lommel Seeliger function 1n our case), and 
then solving for the brightness to assign to a given array element 
location based on the local slope magnitude and direction, relative 
to an assumed sola'^ azimuth and elevation (Batson et al., 1975), 

The gravity data were reduced to Free air and to Bouguer 

r 

anomalies courtesy of the Geopositional Department, Defense Mapping 
Agency Center, St. Louis, Missouri. The reference field at sea • 
level used in the reduction was based on the 1967 International 

Gravity Formula, with all the data referenced to the 1971 

♦ • 

International Gravity Standardization Net. Free air anomalies v/ere 

computed from the following formulation: = g + 0,3086h - y 

Uhere: « Free air anomaly (mgals), g « reading, h = elevation, positive 

down to geoid, y - theoretical gravity value. 

A second order term was added to the elevation correction when the 

magnitude of the correction exceeded 0.1 mgals. Bouguer anomalies 

were computed based oh a slab with density 2.67 gm/cm3 using: 

• ^9BA ** ^9PA ” 0*1119h 
where = Bouguer anomaly in mgals. 

Local terrain corrections were hot included in the Bouguer anomaly- 
computations. 

Generating images from the gravity data presents a more 
complicated problem than the topographic data, since the gravity 
stations are not located along a regular grid. For both the Free 
air and Bouguer anomalies, the data were first scaled to fit within 
a byte variable (i.e., 256 discrete values). The scaled data were 
placed into an array with the same element spacing (30 seconds) and 
geographic coverage as the topography data arr-jy. The element 
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location closest to a given gravity station '' Htion was 

assigned the station value, thereby produci'^ ^ an array occupied in 

part by real data and in part by blank zones. 

We chose to use a conceptually simple, but powerful filtering 

approach for interpolating between the real gravity data points, 

The technique was developed for processing Lunar Consortium data and 

was also used to produce altimetry maps from the aproximately 

100,000 elevation estimates obtained by the Pioneer Venus radar 

« 

mapper (Eliason and SoderOlom,. 1977; Pettengill et al., 1980). The 
technique involves use of a spatial filter that generates a moving 
average through the data array. The filter at any given position 
within the data array occupies a box of N x N elements. The mean 
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■is used as the interpolated value for the box center, only if an 
anomaly value does not already exist at that location. Shifting the 
filter over the array and repeating the operation for each element 
location results in a partially interpolated data set. Blank areas 
still remain in zones without any data and in zones that had so 
little data that a user-defined statistical threshold was not met. 
The threshold is defined in terms of the fraction of elements within 
the box that must have valid data in order to replace the midpoint 

element with the average, assuming that the midpoint did not already 

* » 
have valid data. ‘ 

In practice, the gravity data were processed using several 
filter widths, beginning with a 3 x 3 element filter and ending with 
a El X 21 element filter. The threshold was set such that at least 
20% of the element locations at any given filter position must have 
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had real data values to be interpolated. The choice of filter sizes 

was governed by the station spacings> which averaged several 

kilometers, but varied from 100' s of meters to several 10' s of 

kilometers. A 3 x 3 element filter, for instance, covers 90 seconds 

of arc in both latitude and longitude, or about 3 km along a line of 

latitude. Note that regions with relatively small distances between 

gravity values were interpolated only over relatively short 

distances, since once an element was assigned a data value, it was 

not affected by later filter passes. The net effect of the 

filtering operations was to generate a regular array of gravity 

» • * 

data, with original data left intact, and with interpolated data 

between the original values. Some areas in the array, even with a 

» 

21 X 21 element filter, still had too little data to allow 
interpolation. Those areas appear black in all of cur displays. 

The filtered gravity anomalies could be displayed as gray tone 
or as color-coded images. As with topography, we find that a very 
effective visual display is a shaded relief map. Figures 2 and 3 
depict shaded relief maps for the Free air and Bouguer anomalies, 
respectively. Comparison of these products with previously 
published contour maps (e.g., Woolard and Ooesting, 1964; McGinnis’ 
et al., 1979) show similar broad-scale patterns. However, the maps 
shown in Figures 2 and 3 have inherently more information. The 
reasons are threefold. First, a finer grid spacing was used to 
interpolate between station readings as compared to any published 
map that we have been able to examine. Secondly, there are as many 
contour intervals as there are values in a byte variable (i.e. 256), 


Most contour maps typically have a dozen or so contour intervals 
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Thirdly, areas with a large number of closely spaced stations retain 

details of the anomaly patterns. 

The filtered gravity anomalies can also be digitally overlayed 

onto the topography for a visual display of the correlations between 

the two variables. The technique that we used was to first convert 

the byte value for each element in the gravity arrays to a given 

» 

color (Figure 4). The colors were chosen in such a manner that blue 

• 

corresponded to the most negative anomalies and red to the most 
positive anomalies. Oranges, yellows, and greens were chosen to 
represent intermediate values. Color coded images depicting the 
anomaly patterns were digitally overlayed onto a shaded relief map 
of the topography simply by replacing the brightness component for 
each color with the bniahtness from the corresoonding location in an 
array containing the shaded relief map. The net result is a color 
map, where the brightness is modulated by the shaded relief map 
while the color hue (dominant wavelength) and saturation (degree of 
purity) are controlled by the gravity anomaly values. Note that any 
color can be uniquely defined by a hue (dominant wavelength), 
saturation (purity), and brightness (Slater, 1975). Such products 
are showri for both the filtered Free air and Bouguer anomalies on 
the cover of this issue of EOS. In both cases the maps were 
transformed to a Mercator projection, digitally reduced by a factor 

of seven, and then processed as discussed above. 

» 

DISCUSSION 


A sketch map showing some of .the maj'or structures evident in 
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the color and shaded relief maps Is shown In Figure 5. For this 

* 

paper »ve will restrict our discussion of the processed data to a 

• • 

portion of the midcontinent, where we have a continuing Interest In 

the crustal structure and the relationships between basement 

fractures and ore deposits. Such a discussion also serves to 

Illustrate the kinds of Information that can be gleaned from 

topographic and gravity data that have been processed in image 

» 

format. . . ^ 

The dashed lines in Figure 5 define two major structures in the 
midcontinent, a feature that we call the Missouri gravity low, and 
the Mississippi Valley graben as defined by Kane et al. (1981) on 
the basis of gravity and aeromagnetic anomalies. The Missouri 
gravity low can be seen on the c??lor and shaded relief maps as a HO 
km wide feature that begins at a break near the southern edge of the 
midcontinent gravity high, extends across Missouri, and into the 
Mississippi embayment. In the embayment the low takes the form of a 
horst block called the Pascola Arch (Ervin and McGinnis, 1975). The 
low has an Bouguer amplitude of about 30 mgals for the region just 
to the northwest of the embayment. The magnitude of the anomaly is 
suggestive of a deep basement structure. In fact, the anomaly is 
consistent with an increased crustal thickness of 3 km at the Moho 
or with a crust that has a density of 0.1 gm/cm3 less than the 
surrounding areas’ for the first 4 to 8 km below the surface 
(Strebeck, 1982). Analyses of remote sensing data and ground 
studies for areas overlying the gravity low demonstrate that the 
regional fracture trends are coincident with the strike of this deep 
basement feature, suggesting a fair degree of control over the 
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pattern of faulting in the area (Guinness et al., 1982). 

« 

The intersection of the Missouri gravity low and the 
Mississippi valley graben has clearly served to control the 
emplacement of the two mafic plutons defined by Hlldenbrand ct 
al. (1977): the Bloomington intrusion to thit north and the 
Covington intrusion to the south. The two plutons can be sp/en as 
positive anomalies in the color Free air map and as bumps on ihe 
' shaded relief maps. Additionally, Figure 6 is a gray tone image of 
the filtered Free air data that has been contrast enhanced to show 
the variations in the midcontinent. In such a display, the most 

I 

negative anomalies are dark and less negative anomalies are bright. 

The Missouri gravity low is clearly discernable and there Is a 

suggestion of a NE-Sw trending structufe that, 1n fact, outlines the 

Mississippi Valley graben. The white dots represent about one • 

thousand Earthquake epicenters recorded by the St. Louis University 

• 

seismic network between 1974 and 1979 (Stauder et al., 1977). 

Clearly, the epicenters, which delineate the New Madrid seismic 

high, are concentrated in the crustal block defined by the 

» 

intersection of the Missouri gravity low and the Mississippi Valley 

• • 

graben. Thus, the intersection of these two crustal structures 
seems to have provided both a zone of weakness for the emplacement 
of plutons, and a focus for release of strain energy via 
Earthquakes. 

Examination of the Free air color map suggests an extension of 
the Missouri gravity low from the midcontinent gravity high to the 
northwest, perhaps as far as to the Big Horn uplift in Wyoming. The 
suggestion is based on an alignment of a sharp break in the Free air 
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anomalies with the North Platte River valley in Nebraska^ and with 
fractures in Wyoming. The reality of such an extension must remain 
open to question until further study. If the extension is 
genetically related to the Missouri gravity low, the combined 
lengths of the two features would be over 1500 km. There are other 
structures of comparable length on the Earth, namely transcurrent 

faults such as the Altyn Tagh of China (Molnar and Tapponier, 1975). 

* 

It is conceivable that the Missouri gravity low and the extension 
into Wyoming are part of a transcurrent fault system, sections of 

ft * 

which have been reactivated during various time periods. Such an 
origin is also consistent with the observation that the structure 
cuts across the basement age and province boundaries defined by 
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for instance, during the postulated collision event that produced 
the Grenville orogeny (Dewey and Burke,’ 1973). 


IMPLICATIONS FOR THE EARTH SCIENCES 

The future will see an increased availability of a number of 

* 

other large, geographically-oriented data sets that are applicable 
to Earth Science problems. For instance, as part of the National 
Uranium Resource Evaluation Program (NURE), the Department of Energy 
has acquired airborne digital gamma ray (net radioactivity) data, 

t 

aeromagnetic anomalies, and hydrogeochemical data for most of the 
continental United States. Such large data sets are aptly suited 
for use with the kinds of techniques described in this paper. Of 
course, analysis of any data set requires above all the proper, N 
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framing of scientific questions. Also* the right kinds of data must 
be accessible to answer the questions. However, we do feel that 
interactive digital image processing of a variety of 
geographically-oriented data sets can lead to substantial advances 
In our understanding of the Earth. 
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Figure 1 


Figure 2 
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FIGURE CAPTIONS 


- Shaded relief map depicting NOAA 30 second elevation 

t > 

averages for the continental United States, Numbers 

running vertically along the sides of the image are 

degrees of north latitude, while numbers along the bottom 

are degrees of v/est longitude. This and other maps in 

the paper are Mercator projections. The simulated sun 

west was set at 20* above the western horizon. Dark • 

blocky areas are missing data. 

* 

~ Shaded relief map depicting filtered Free air anomalies, 
Black areas indicate regions containing either bad data 
or too little data to have a statistically valid 
interpolations. The simulated sun was set at 15* above 
the western horizon. Variations in the anomalies are, in 
effect, illuminated as if they were a set of topographic 
ridges and valleys. For example* the midcotinent gravity 
high (43*N.» 94*W.), which is a positive anomaly, stands 
as a ridge, while the great valley of California (37*N., 
120*W.), which is a negative anomaly, appears as a 
trough. 

- Shaded relief map depicting filtered Bouguer anomalies! 
The map was produced in the same manner as the Free air 

. anomaly map in Figure 2. No local terrain corrections 
were, included in computing Bouguer anoma’iies. 
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Figure ^ - This diagram shows the approximate color-coding scheme 

that was used to produce the color contours for the 

gravity anomaly maps shown on the cover of this issue of 

EOS. Histograms showing the fractional area occupied by 

each anomaly value were used to assign colors in such a 
» • 

manner that each color value occupies an area comparable 
to any other color value. About 45 discrete color values 
(i.e,, discrete hue and saturation values) were used in 

I 

each anomaly map. On the other hand, the shaded relief 
map versions of the anomalies each contain ?56 shades of 
g»'ay. 


Figure 5 - Sketch map showing some of the major structures that can 
be seen in the gravity anomaly maps presented in this 
paper. The Missouri gravity low can be seen best' on the 
shaded relief map versions of the anomalies. The two 
plutons are the Bloomfield pluton to the north and the 
Covington pluton to the south. They can be seen as red 
dots on the Free air color maps and as bumps in the 
shaded relief anomaly maps. The Bloomfield intrusion is 
located at about 37*N., 90*W. The Mississippi valley 
qraben location is from Kane et al. (1981), although the 
* outline can be seen in the shaded relief anomaly maps. 

Figure 6 - Gray tone version of the filtered Free air anomah’es. 

The map boundaries extend from 34* to 43*N.» 87*. to 99*W. 

The two horizontal lines correspond to 37* and 

40*N. latitude, while the two vertical lines are 91* and 
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95*W. longitude. In this map bright areas have large 
Free air anomalies while dark areas have small (i.e., 
highly negative) values. The white dots are Earthquake 
epicenters for the period from 1974-1979 as recorded by 
the St. Louis University seismic network. The Bloomfield 
and Covington plutons are the bright circular features 
, just to the north and south, respectively, of the 
epicenter cluster. The Missouri gravity' low extends from 
a break in the midcontinent gravity high, through 
Missouri, and across the Mississippi valley graben. The 
intersections of the gravity low and the graben are the 
.locations of the intrusions. In addition, most of the 
seismicity is concentrated within the crustal block • 

defined by the intersection of the two gravity low and 

* , 

the graben, 


CAPTION FOR COVER:* 


About half a million Free air and Bouguer anomalies for land 

stations covering the United States were spatially filtered to form 

a regular array of interpolated data. The anomaly values were then 

color-coded and overlayed onto a shaded relief map of topography. 

» 

The simulated sun for the shaded relief map is from the west at 20“ 
above the horizon. The color-coding scheme for the Free air data is 
approximately as follows: blue, less than -34 mgals; blue-green, 
-34 to -10 mgals; green -10 "to -3 mgals; yellow, -3 to +10 mgals; 
orange, +10 to +26 mgals; red, greater than +26 mgals. The 
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color-coding scheme for the Bouguer data is approximately as 
follows: blue, less than -212 mgals; blue-green, -212 to -128 
mgals; green -128 to -84 mgalsj yellow, -84 to -47 mgals; orange, 
-47 to -17; red, greater than -17 mgals. See the article in this 
issue of EOS for further information. 
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